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Abstract 

Coherent anti-Stokes Raman scattering (CARS) microspectroscopy has demonstrated significant potential for 
biological and materials imaging. To date, however, the primary mechanism of disseminating CARS spectroscopic 
information is through pseudocolor imagery, which explicitly neglects a vast majority of the hyperspectral data. 
Furthermore, current paradigms in CARS spectral processing do not lend themselves to quantitative sample-to- 
sample comparability. The primary limitation stems from the need to accurately measure the so-called nonresonant 
background (NRB) that is used to extract the chemically-sensitive Raman information from the raw spectra. 
Measurement of the NRB on a pixel-by-pixel basis is a nontrivial task; thus, reference NRB from glass or water are 
typically utilized, resulting in error between the actual and estimated amplitude and phase. In this manuscript, 
we present a new methodology for extracting the Raman spectral features that significantly suppresses these 
errors through phase detrending and scaling. Classic methods of error-correction, such as baseline detrending, are 
demonstrated to be inaccurate and to simply mask the underlying errors. The theoretical justification is presented 
by re-developing the theory of phase retrieval via the Kramers-Kronig relation, and we demonstrate that these 
results are also applicable to maximum entropy method-based phase retrieval. This new error-correction approach 
is experimentally applied to glycerol spectra and tissue images, demonstrating marked consistency between spectra 
obtained using different NRB estimates, and between spectra obtained on different instruments. Additionally, in 
order to facilitate implementation of these approaches, we have made many of the tools described herein available 
free for download. 


Introduction 

Coherent anti-Stokes Raman scattering (CARS) is a nonlinear scattering phenomenon in which two photons, “pump” 
and “Stokes”, coherently excite a molecular vibration. From the excited mode, a “probe” photon inelastically 
scatters off with an energy increase equal to that of the vibrational state. This optical technique affords label-free, 
molecularly-sensitive investigation of samples without autofluorescence competition and at significantly higher speeds 
than offered by traditional spontaneous Raman spectroscopy[l, 2]. CARS microscopies and microspectroscopies have 
demonstrated success on a large variety of material and biological samples ranging from polymer blends[3, 4] to brain 
tumor masses [5, 6], requiring fractions of a second for a micrograph highlighting a single vibrational mode to a few 
minutes for a complete hyperspectral image. 

A fundamental challenge to harnessing the information content in CARS microspectroscopies, is extraction of 
the chemically-specific Raman signal from the so-called “nonresonant background” (NRB). The NRB is predomi¬ 
nantly comprised of electronic signal contributions from other nonlinear optical phenomena that are less chemically 
specific)?]. Although it is sometimes viewed as an interference, the NRB actually amplifies the weak Raman sig¬ 
nal, enabling high-sensitivity detection]?, 6]. This coherent mixing, however, does bring spectral distortions. For 
CARS microscopies that probe small increments of the full vibrational spectrum, physical methods are utilized to 
reduce the NRB generation, which in turn reduces the overall CARS signal [8, 9, 10, ?]. For spectroscopic CARS 
techniques (microspectroscopies), however, two classes of numerical methods are commonly used to remove the dis¬ 
tortion of the NRB: one based on maximizing entropy [11] and the other utilizing the Kramers-Kronig (KK) relation 
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[12]. These techniques are functionally equivalent [13] and have demonstrated success with materials[14, 15] and 
biological samples[16, 17, 18, 19, 6]. Both of these techniques, however, rely upon an accurate measurement of the 
NRB, which also incorporates the system excitation profile and spectral response. To date, no approach has been 
found to conveniently and accurately measure the NRB, so retrieved spectra generally have phase and amplitude 
errors. Additionally as typically implemented, these techniques implicitly assume a constancy of the nonresonant 
background at every pixel across the image, which is typically a poor assumption for heterogeneous samples. The 
resulting processed spectra contain baseline fluctuations, distorted peaks, and provide limited quantifiable informa¬ 
tion between molecular species; thus, intra- and inter-sample comparisons are significantly hampered. Additionally, 
the measurements performed on different CARS instruments- or even the same instrument under altered operating 
conditions- are not directly comparable. 

In this work, we present a new method of processing CARS spectra that suppresses errors resulting from the use 
of an inexact reference NRB spectra, removing baseline fluctuations and generating spectra that are agnostic to the 
reference material used. In upcoming sections, we present a theoretical re-examination of the KK relation and errors 
resultant from reference NRB measurements. Use of the KK presents analytical expressions for correcting these 
errors, but these results are also broadly applicable to the maximum entropy method (MEM)[II]. Finally, we present 
simulated and experimental results using neat glyercol and murine pancreas tissues. Spectra are presented from two 
broadband CARS (BCARS) architectures with dramatically different excitation profiles[I9, 20, 6]. The underlying 
purpose of this manuscript is to enable the extraction of “pre-processed” spectra that are universally comparable in 
amplitude and shape. This facilitates dissemination of pre-processed spectra as a new “currency” of coherent Raman 
imaging data. 


1 Theory 


Classically, the CARS spectral intensity, Icars^ may be described as[21]: 


IcARs{^) oc 


{uj)Ep{ujp)Eg{uJs)Epr{ujpr) X S {uj - UJp + UJs - ujpr)dujpdujsdujpr) 


( 1 ) 


where is the nonlinear susceptibility of the sample; Ep, Es, and Epr are the pump, Stokes, and probe field 
amplitudes, within the frequency spaces, tOp, lor, and ujpr, respectively; and the delta function ensures energy 
conservation. This equation may be re-written in a more intelligible form[6]: 


IcARsi^) OC 


[Esiuj) A Ep{uj)] x^^\uj) * Epriuj 

Cst 


( 2 ) 


« (3) 

where * and * are the cross-correlation and convolution operations, respectively, and Cst is the coherent stimulation 
prohle. Equations 1 and 2 are mathematically equivalent. If we assume a spectrally narrow probe source, we 
can introduce an “effective” stimulation profile, Cst, and nonlinear susceptibility, as presented in eq 3, where 
Cst (to) = [Cstioj) * Epr{(jj)] / / Epr{uj)dw and * Epr (to). 

The nonlinear susceptibility describes signal contributions from Raman vibrationally-resonant, XR, and vibra- 
tionally non-resonant, xnr, sources: + XNRi^^)- To a hrst degree approximation[22], spon¬ 

taneous Raman spectra, IRaman, are related to the vibrationally resonant component of the CARS spectra as: 
lRaman{<^) OC Im{xR(w)}, where ‘Im’ indicates the imaginary component. The purpose of phase retrieval is to 
ascertain a phase, that isolates the Raman resonant components from the total nonlinear susceptibility. 


1.1 Phase Retrieval Using the Kramers-Kronig Relation 

The KK relation states that there is an explicit relationship between the real and imaginary components of a function, 
/(w), that is causal (analytic[23, 24]); thus, if only the real (or imaginary) part is known, the imaginary (or real) part 
can be calculated. In CARS and other spectroscopies, neither the real nor imaginary portion of x^^^ is accessible {n.b.: 
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Cst in eq 3 is not a causal function). If the function is square-integrable, there also exists an explicit relationship 
between the complex norm of the function and the phase[25]: 

ln|/(w)| = -n{(j}iuj)} (4) 

(j){uj) = 'H{ln|/(w)|}, (5) 


where Ji is the Hilbert transform[26], which assumes knowledge of the complex modulus of the function or phase 
over an infinite frequency range. Practically, however, the CARS spectral recording window is limited; thus, we will 
introduce a “windowed” Hilbert transform: TLw, as: 

nw{j(x);u}a,i^b} = — / - jdx 

^ JuJa ^ 

lim 'Hw{fix);u;a,ujb} = 'H{fix)}, 

UJa—^ — OO 

CJb—^OO 

which is limited to the spectral range Wa to ujb (for compactness, these parameters will be neglected from the operator 
form). V is the Cauchy principle value. Under the conditions that (a) the Raman peaks encompassed within this 
window are not affected by those outside of the window and (b) any resonances of xnr are far-removed from those 
of xr (as is typically the case with infrared stimulation), the windowed- and analytic Hilbert transform are related 
as: 

|^ln|x^^Hw)p| « ?i|iln|x(3)(a;)|2| (8) 
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where e{uj) is an additive error term (see Figure SI and associated Supporting Information demonstrating the additive 
nature). Applying eqs 5 and 8 to eq 3, the retrieved phase from the raw CARS spectrum, 4>cars, niay be described 
as: 
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where Z denotes the angle (phase). The retrieved phase is not simply that of the nonlinear susceptibility but 
also contains contributions from the windowing error and the effective stimulation profile. If one can measure the 
NRB spectrum without Raman components, Inrb, and assuming that the spectrum is far-removed from electronic 
resonances such that xnr is approximately real, then the following phase retrieval can be utilized in lieu of eq 10: 

(t>CARS/NRB{^) = Tdw 12 InrbIm^) } ^ | 2 } 

— e(a;)-b "Hiv In |C'st(w)p| 

+ ^ [Xfl(w) + XJVfl(w)] - ZxNRiuj) 

« ^[Xfl(w)+XJVfl(w)], (11) 


which is analogous to applying the KK relation to Icars/Inrb- Using this ratio as our signal, the retrieved complex 
spectra, Iretr, is: 


„(w) = 


' Icars{< jj) 
Inrb{^) 


exp i(j)CARS/NRB 
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and the Raman-like spectrum, Irl: 


Irl{<jJ) =Im{Aeir(w)} 


Ini{Xfl(^)} . 

IXAffll 


(13) 


thus, the Raman-like spectrum is proportional to the spontaneous Raman spectrum scaled by the nonresonant com¬ 
ponent. The earliest KK relation results [12] (presenting a very different derivation) utilized VIcars saa'Pcars/nrB) 
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which is directly proportional to the spontaneous Raman spectrum, but it implicitly assumes that Cst (w) is constant 
and does not account for e(a;). Later work[18], using an identical phase retrieval method (although with a different 
derivation- see the Supporting Information), identified a need to remove the stimulation profile via Icars/Inrb- 
Together these works explain the method of phase retrieval under certain ideal conditions. In the following sections, 
we will present the ramifications of the real-world condition when the NRB of the sample is not directly measurable. 
This analysis, as only enabled by the derivations in eqs 8 through 13, will lead to a direct method to account and 
correct for these errors. 


1.2 Errors from Inaccurate NRB Measurement 


Although phase-retrieval is theoretically straightforward, measuring the NRB is technically challenging. To circum¬ 
vent this limitation, researchers rely on measurements from model materials, such as coverslip glass or water. This 
practice assumes that the model material does not contain a significant Raman signature and the NRB varies little 
from material to material. Additionally in practice, any reference material spectral deviations from the actual NRB 
are assumed to contribute to a slowly-varying baseline that can be subtracted off via various detrending methods. 
With the increased sensitivity of spectroscopic CARS systems, however, these assumptions appear more and more 
naive. As demonstrated below, the difference between the NRB and reference measurement does not lead to an 
additive error, but rather a multiplicative complex error. 

We obtain a reference measurement, Iref, as a surrogate for a proper NRB measurement. Here, Jre/(w) = 
5(w)/Ar_R_B(w), and ^(w) is assumed to be real and positive. By applying eq 11: 

(14) 


4 ’err (^) 


4’CARS/ref{‘^) — <PcARS/NRb{a>) +'Hw 


-In^ 

2 C(w) 


we see that the Raman-like spectrum (eq 13), I^L-ref, is: 

^RL—ref{^') — 
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From eqs 14 and 15, one can see that the use of a reference measurement leads to both amplitude (Aerr) and 
phase {(j)err) distortions. Accordingly, removing these errors is not simply a matter of subtraction. The phase error, 
however, is additive in nature, and connected to the amplitude error via the KK relation: 


^err(a^) — ^ {la Ag,.,. (cc ) } 
In A 

err (w) = —'H{4> err (w)}- 


(16) 

(17) 


There is, however, an ambiguity in this relationship. If ^(w) is multiplied by a constant, (perr {uj) = 'H{lnl/.=,^(w)} = 

77 {lnl/^(w)}, since the Hilbert transform of a constant is zero. 


1.3 Correcting Phase Error and Scale 

The purpose of correcting for Raman signature extraction errors is not simply to generate qualitatively accurate 
spectra, but those that are quantitatively reliable, facilitating direct comparison and analysis of spectra collected 
of different samples with potentially different reference materials, and on various spectroscopic systems having dif¬ 
ferent excitation profiles. In the previous subsection we demonstrated that the use of a reference NRB which only 
approximates the nonresonant response of the material induces amplitude and phase distortions. Additionally, the 
commonly used method of subtracting baseline fluctuations, borrowed from the spontaneous Raman community, 
does not actually remove the errors as (1) the nature of the error is complex valued and (2) the amplitude error is 
multiplicative. 

Below we show that one can properly correct for signal extraction error using the following steps: 

1. Remove phase error via phase detrending, and correct for part of the amplitude error via the KK relation 

2. Correct for scaling error (involving S) and the error from the windowed Hilbert transform (of (perr) via unity¬ 
centering of the real component of the retrieved (phase-corrected) spectrum. 
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As displayed in eq 14, the difference between the ideal phase retrieval (in which the NRB of the sample is exactly 
known) and that using a model material is (f)em which is additive. The retrieved phase (ideal) is qualitative similar 
to Raman-like spectra in that the spectral features are peaks that extend positively from a zero baseline. A slowly- 
varying phase error will cause a slowly-varying deviation from the zero baseline. Finding (j)err', therefore, is a matter 
of isolating the erroneous baseline, for which one may use traditional baseline detrending methods (albeit, applied 
to the phase rather than the Raman-like spectral amplitude). From this extracted (j)err^ using eq 17, one can find 
the amplitude error. With these variables in hand, one can multiply the retrieved complex spectrum by a complex 
phase-correction multiplier, generating a phase-corrected spectrum, Ipc'. 


I pc (^ ) - 


' IcARsi<^) 

/re/(w) 


exp [i4>CARS/refi‘^) 


exp -y. {4>erril^)} 


exp [-t(()err(w)] 


(18) 


n.b.: the use of “phase-corrected” in this context is unrelated to that in Ref. 18. As described in the previous 
subsection, however, there is an ambiguity in the KK relation between the error phase and amplitude. Thus the 
calculated A^rr is only accurate to within a constant multiplier. Additionally, the Hilbert transform in eq 18 is 
actually a windowed Hilbert transform; thus, H {(fierri^jj)} = T~Lw {</>err(w)} -I- eerr(w), where Cerr is a window-effect 
error term similar to that introduced in eq 8. 

To finalize the error correction, one needs to account for the A^rr ambignity and terr- Both of these quantities 
are discoverable by analyzing the real component of the phase-corrected spectrum in eq 18 since the real component 
of eq 12 is unity-centered, i.e., (Ix^^^I/Ixat^rI cos4>CARs/nrb) = 1- The existence of S, however, will alter the mean; 
thus, one could calculate the mean of the real component of the retrieved spectrum and normalize the real and 
imaginary components by this value, tem however, may impart a frequency-dependent component to this mean. 
Using numerical means, though, one can find a slowly-varying center-line and normalize the phase-corrected spectrum; 
thus, removing S and terr in one step. Assuming this centerline can be found, a complex corrected spectrum Icorrected 
may be calculated: 


^corrected 


(w) = 


Tpc (^) 


(Re{/pc(w)}) (w) 

-j, ^ ^ 

—-exp 1(PCARS/NRB (w), 

\Xnr{uj)\ 


(19) 

( 20 ) 


where we have noted the frequency-dependence of the mean-line in the denominator of eq 19. Comparison of eqs 12 
and 20 shows that using the prescribed steps, one can retrieve the same spectrum using a reference NRB as if the 
NRB were measurable. 


2 Materials and Methods 

2.1 CARS Microspectroscopy 

In this manuscript, most CARS spectra were collected on a recently introduced broadband CARS (BCARS) system 
that has been described elsewhere [6]. In order to demonstrate that properly retrieved spectra can be essentially 
identical, irrespective of instrumentation, some spectra were collected on an instrument that uses an earlier im¬ 
plementation of BCARS signal generation [19, 14]. The newer BCARS system excites molecular vibrations more 
efficiently with the highest response at the lowest wavenumbers, whereas, the traditional BCARS system excites 
most Raman transitions with relatively uniform response. For the newer BCARS system, the total average incident 
power was <25 mW (3.5 ms integration time), and for the traditional BCARS system <60 mW (7.8 ms integration 
time). 


2.2 CARS Simulations 


The CARS simulation software was developed in MATLAB (Mathworks) and numerically implements eq 2 directly. 
All excitation sources were simulated as real Gaussian functions, and the Raman response a complex Lorentzian 
(damped harmonic oscillator) as: 


xfl(w) = 

m 


VLr. 


Am 

-UJ - iTm ’ 


( 21 ) 
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where Am, ^m, and Tm describe the amplitude [multiplier], wavenumber, and half-width of the Raman peak. 
Specific parameters for the simulations are provided below. 

2.3 Signal Pre-Processing 

Image and spectral processing are performed using in-house software developed in MATLAB. Specific details of the 
pre-processing of BCARS spectra (or images) are presented in the Supporting Information (see Figures S2 - SII 
and Table SI). In brief, dark spectra are collected as are NRB spectra from reference materials (glass coverslip, 
glass microscope slide, or water). For noise reduction of BCARS hyperspectral data, we utilize singular value 
decomposition (SVD)[27] on Anscombe transformed spectra. The Anscombe transform[28, 29] normalizes the noise 
variance, accounting for mixed Poisson-Gaussian noise. Pertinent singular values are selected by noise analysis in the 
spectral and spatial domains in an automated or semi-automated fashion (Figures S8, S9). Once SVD is performed, 
the variance-stabilized, noise-reduced spectra are returned to their normal mixed-noise state using an optimized, 
generalized inverse Anscombe transformation[29]. 

Raman-like spectra are retrieved using the Hilbert transform implementation of the KK relation (described later). 
The erroneous component of the retrieved phase is found in an automated fashion using an asymmetric lease-squares 
(ALS) technique with a Whittaker smoother[30, 31, 32]. Phase and partial amplitude correction are performed as 
described in eq 18. To determine the mean trend line for final spectral correction (eq 19), a Savitky-Golay filter is 
utilized. 

2.3.1 Phase Retrieval Using the Hilbert Transform 

The Hilbert transform (eq 6) is implemented in the time-domain (t): 

Hw {/(w)} =T{i sgn(t) {/(cc)}} , (22) 

where T and are the Fourier and inverse Fourier transforms, respectively, sgn(t) is the signum (“sign”) function 
and /(w) is a spectrally dependent function (e.g., Icars / Iref)- Additionally, we implement a spectral padding 
procedure [11] to extend the window range, reducing numerical edge effects. This method efficiently retrieves phase 
with only two Fourier transforms and was designed for parallel processing. 100 parallel solutions, for example, 
with each spectrum containing 1000 spectral points requires ~200 /rs per spectrum on a personal computer (16 
GB RAM, 3.4 GHz quad-core processor). The KK and associated Hilbert transform code is freely available at 
http://github.com/GoherentRamanNIST in the MATLAB and Python languages. 

3 Results 

3.1 Simulated Spectra 

To validate the presented theory on phase retrieval and error correction, we begin with the simplified case of a 
two-peak Raman system with parameters (eq 21): Ai = 0.25, Hi = 1000 cm“^, Fi = 10 cm“^, A 2 = 1, H 2 = 3100 
cm~^, and F 2 = 20 cm”^). xnr = 0.55 and Xref is Xnr multiplied by a Gaussian function. The simulated nonlinear 
susceptibilities are presented in Figure 1 a. The GARS spectra that result from these susceptibilities are shown in 
Figure 1 b. 

Figure 2 a shows the retrieved spectra using a reference or the actual NRB (“ideal”). The reference-retrieved 
spectrum show a clear, large baseline and distorted (asymmetric) peak shapes. Additionally, the peak amplitudes are 
sufficiently perturbed that the 1000 cm“^ peak appears ~50% larger than the 3100 cm“^, although the latter should 
be the larger of the two. Figure 2 b shows the difference (A) between the ideal and nonideal retrieved spectra. From 
this, one can gather that the traditional tactic of baseline detrending will only resolve the slowly-varying baseline, but 
the underlying peak-errors (amplitude and phase) will remain. Figure 2 c shows the phase retrieved under ideal and 
nonideal conditions, which does not display any obvious peak distortions. As clearly presented from the difference 
of the retrieved phases (Figure 2 d), there is no spectral distortion of the Raman peaks, but only the slowly-varying 
baseline {4>err in eq 14). 

Using this phase error and applying a calculated amplitude correction (eq 18), the baseline and asymmetric 
spectral distortions are removed entirely, as shown in Figure 3 a (differences shown in Figure 3 b). For reference, a 
traditional amplitude detrending is also displayed, showing the remaining distortions clearly. Although the phase- 
corrected spectrum is qualitatively similar to the ideal, the relative amplitudes of the two peaks are still incorrect, 
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Figure 1: Simulated nonlinear susceptibilities and CARS spectra, (a) Raman {xr) nonresonant nonlinear 
susceptibilities of the sample (xnr) and of a reference {Xref)- (b) Simulated CARS spectra of the sample {Icars, 
X = XR + Xnr), the sample in absence of Raman components {Inrb), and of the reference material (Iref)- 


owing to Cerr and the amplitude ambiguity described previously. Figure 3 c shows the real part of the phase-corrected 
spectrum and the mean trend-line that deviates from unity. Using this trend as a scaling factor (eq 19), Figure 3 d 
demonstrates that the phase-corrected spectrum is now identical to the ideal retrieval (difference shown in Figure 3 
e). 

As presented elsewhere[13], the KK and MEM phase retrieval methods are functionally equivalent. Figure S12 
demonstrates the applicability of the presented phase error correction method to Raman-like spectra extracted from 
simulated BCARS spectra via the MEM method. Like the KK demonstration above, the prescribed method enables 
a reference NRB spectrum to be utilized and generates a corrected Raman-like spectrum that is equivalent to one 
extracted when the NRB is exactly known. 

3.2 Experimentally Measured Spectra 

The developed error correction method can readily be applied to experimental results without modification. Ad¬ 
ditionally, this method provides spectra that are comparable between microscopy platforms. Figure S13 a shows 
CARS spectra collected for neat glycerol on two BCARS systems (average of 100 acquisitions), demonstrating widely 
different system responses for the same molecule. Figure S13 b shows the recorded reference spectra for 3 different 
commonly used model materials. These reference spectra not only demonstrate a great variety of overall amplitude 
but also spectral features. As expected, retrieving the Raman-like spectra using these references produces amplitude 
and phase errors, resulting in distortions as shown in Figure S13 c. Figure 4 a shows the Raman-like spectra with 
the slowly-varying baseline removed. The spectrum retrieved using water demonstrates the most severe distortions. 
Even the spectra using glasses (slide and coverslip) demonstrate differences. In comparison. Figure 4 b shows the 
same four spectra after full correction (eq 19). The four spectra are significantly more similar in amplitude and 
shape. Additionally, one should notice the partial recovery of the OH-stretch peaks (^3300 cm“^) for the spectrum 
retrieved using water. When a particular reference material is utilized, the retrieved spectra will have suppressed 
peaks wherever the model material contains Raman peaks. Within spectroscopic CARS literature, coverslip or slide 
glass have often served as a convenient reference. What was not apparent at the time, however, is that these glasses 
have nontrivial, glass-dependent Raman peaks. The primary cause of the spectral differences in Figure 4 b are 
due to these reference material Raman peaks. Figure S15 b shows the retrieved (and corrected) spectra of these 
reference materials, with their peaks exactly correlating with the deviations in the retrieved spectra (Figure S15 
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Figure 2: Comparison of retrieved spectra using the KK under ideal conditions and with use of am NRB reference, 
(a) Retrieved spectrum (Iretr) when the NRB is exactly known (Ideal) and when a reference is utilized (Ref), (b) 
The difference, A, between the ideal and reference retrieval showing remnants of the original Raman peaks, (c) The 
KK-retrieved phase when the NRB is known {4>CARs/nrb) and when using a reference {4>CARS/ref)- The phase 
difference, A, between the ideal and reference retrieval that shows only a smooth baseline with no Raman peak 
information. 


a). These spectra were collected using a time-windowing, self-referencing method (see Supporting Information and 
Figure S14). This enables collecting an NRB-dominant spectrum directly from the sample, which can then be used 
as the reference. This also enables the use of reference material spectra with their Raman peaks suppressed. Figure 
S15 c shows Raman-like (corrected) spectra of glycerol using different references with their peaks suppressed. 

3.3 Tissue Imaging 

The presented phase retrieval method can also be reliably applied to hyperspectral images. For this purpose, we 
imaged a histological section of murine pancreatic artery. A 200 /rm x 200 section (90,000 pixels total) was 
imaged with 3.5 ms dwell times. Reference spectra were collected from water and the sample coverslip. The raw 
BCARS image was de-noised using SVD on Anscombe-stabilized spectra, keeping 23 singular values. After this 
de-noising, the hyperspectral data were processed four times: twice with each reference spectrum and twice with 
amplitude or phase detrending methods. Figures 5 (a,c) are the pseudocolor images of the murine skin highlighting 
protein in blue (2937 cm“^- 2882 cm“^), DNA in orange (785 cm“^), and in red a shoulder peak that is dominant 
in the smooth muscle (1339 cm“^) and is tentatively assigned to actin/myosin[33, 34]. The Supporting Information 
provides more detail regarding how the images were processed. The left-halves of Figures 5 (a,c) show the processing 
performed using the coverslip NRB reference, and the right-halves using water. The color intensity range of red, 
green, and blue channels are the same for both half-images within a single hgure; although, the range is different 
between Figure 5 a and Figure 5 c. With amplitude detrending (Figure 5 a), the boundary between the half-images is 
obvious. Using the coverslip reference, the blue channel (protein) is more intense than when using water. Conversely, 
the red channel (actin/myosin) is suppressed when using coverslip. The DNA/RNA signature is nearly the same 
in both images. This highlights that distortions vary across the spectrum, and one cannot simply normalize out 
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Figure 3: Residual error and their correction using the proposed phase-detrending and scaling with comparison to 
tradition amplitude detrending, (a) KK-retrieved spectra when the NRB is known (Ideal) and using a reference 
material with phase detrending (Ph. Det., see eq 18) and amplitude detrending (Amp. Det.). (b) The difference, A, 
between the ideal and corrected spectra. Amplitude detrending creates asymmetric peak distortion not present with 
phase detrending, (c) Using the mean-trend of the real component of the retrieved spectrum allows for correction 
of edge effects and scaling ambiguity, (d) Phase-detrended and scaled spectrum (Ph. Det. + Scaling, imaginary 
portion of eq 19) is identical to the ideal retrieval (Ideal), with (e) no difference, A. 


the errors. Using phase detrending (and scaling) in Figure 5 b, there is no obvious discontinuity between the two 
half-images. Figure 5b shows single-pixel spectra collected from a portion of the internal elastic lamina (marked by a 
white ’x’ in Figures 5 a) and corrected using amplitude detrending. In these spectra the coverslip-processed spectrum 
is ^40% stronger at higher energies, but ^15% to 30% weaker within the fingerprint region. With phase detrending, 
shown in Figure 5 d, the peak differences are significantly less obvious, with the predominant cause of residual error 
being from Raman peaks inherent to the reference materials (in Figures 5 (b,c) significant perturbations induced 
by reference NRB Raman peaks in coverslip glass are denoted with a and dashed lines). Figure S16 shows a 
histogram comparison of the Raman peak amplitudes whether performing amplitude detrending or phase detrending 
(and scaling), demonstrating significantly increased similarity when using the developed method. To quantify the 
improvement, we calculated on a pixel-by-pixel basis the relative difference (Figure 5 e) between the peaks used to 
construct Figures 5 (a,c), i.e., the coverslip-processed intensity minus the water-processed intensity, over the mean. 
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Figure 4: Spectral correction of neat glycerol, (a) KK-retrieved Raman-like spectra using NRB reference materials 
with traditional amplitude-detrending, highlighting the inability of baseline detrending to remove amplitude and 
phase error, (b) Corrected spectra with phase detrending and scaling (imaginary portion of eq 19), showing close 
agreement with residual differences primarily arising from reference material Raman peaks. 


Using amplitude detrending, the mean relative difference is 10.5% (standard deviation, tr, 3.7%) for protein, -17.4% 
(cr=8.4%) for actin/myosin, and -1.1% {a = 0.5%) for DNA/RNA. Using phase detrending and scaling, the mean 
relative difference for protein improves to 3.5% {a = 4.5%) and to 0.08% (cr = 1.1%) for actin/myosin. The mean 
relative difference for DNA/RNA remains nearly the same (in amplitude) at 1.7% (cr = 0.5%). 


4 Summary and Discussion 

The aim of this work is to foster quantitative reliability and repeatability of CARS spectra and to promote quantitative 
analysis of hyperspectral data cubes rather than rely on the pseudocolor imagery. Although this work presents the first 
efforts towards making CARS spectra reliable, repeatable, and universally comparable, there are still many avenues 
to improve upon this method. For example, the retrieved spectra using this (and other) methods are normalized to 
the NRB. As described in eq 12, this results from the practical necessity of removing the system response, which 
may contain spectral fluctuations due to laser source spectral profiles or hlter set characteristics. One solution may 
be a thoroughly characterized reference material of which the resonant and nonresonant nonlinear susceptibilities 
are known; thus, enabling the extraction of the system response. This extracted system response could then be 
applied to the analyses of future CARS acquisitions on that system. Another opportunity is the examination of the 
proposed method under the condition that the nonresonant nonlinear susceptibility is not approximately real, but 
rather complex. The ramifications of a complex NRB may become more appreciable due to water absorption as 
CARS systems move further into the infrared or due to multiphoton electronic absorption as more exotic fluorescent 
species such as plants or algal cells are investigated. 

In conclusion, this work presents an expanded method of Raman signal extraction that corrects for amplitude 
and phase errors that are ubiquitous in CARS microspectroscopy. The presented technique, does not require any 
augmentation to current acquisition work flows and performs the corrections in silico. The corrected spectra, as 
demonstrated with neat liquids and tissues images, significantly reduces the intra-spectral distortions caused by 
the use of NRB reference spectra, and facilitates direct, quantitative comparison between samples and microscopy 
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Figure 5: Spectral retrieval and correction of hyperspectral data: murine pancreatic artery tissue section, (a) 
Traditional method of preparing pseudocolor imagery shows distinct differences when using a glass coverslip or water 
as the NRB reference, which cannot be simply corrected with normalization, (b) Single-pixel spectra (white ’x’ in 
(a)) show intra-spectral deviations that are reference-dependent. Additional errors due to Raman peaks emanating 
from the reference NRB materials (“*’ and dashed lines), (c) Pseudocolor imagery using the presented phase retrieval 
shows no obvious sign of difference between the halves processed with different NRB references, (d) Single-pixel 
spectra show close agreement with residual error due to Raman peaks emanating from the reference NRB materials 
(‘*’ and dashed lines), (e) Histogram analysis of the relative difference between spectral peak intensities used in the 
creation of (a). The developed method shows <4% difference (mean) between peak amplitudes. 


systems. In the future, it is the hope of the authors that this method (and future improvements) will enable mass 
dissemination of coherent Raman hyperspectral data cubes for community data mining and analysis. 
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SI Theory 

Sl.l Additive Nature of Window-Edge Effects with the Numerical Hilbert Transform 

Within the main text of the manuscript it was proposed that the full Hilbert transform and the “windowed” Hilbert 
transform were related by an additive error term, e, under the condition that the resonant component of the nonlinear 
susceptibility, XRi fully captured and the nonresonant component, xnr, was not. That is: 

{loglXfi(w)+ Xvfl(w)|} = {loglxfl(w)+Xvfl(w)|} + e(w) (SI) 

= ^ [Xfi(w) + Xvfl(w)], (S2) 


where Z is the phase. 

To demonstrate this relationship, we simulated a nonlinear susceptibility with two Raman peaks and an extremely 
broad (25,000 cm“^ full-width half-max) nonresonant term as shown in Figure SI a. The windowed Hilbert transform 
of the complex norm of the nonlinear susceptibility was performed over larger-and-larger windows and compared with 
the ideal solution (the phase). Figure SI b shows the retrieved phase via the windowed Hilbert transform with varying 
window widths. Figure SI c shows the numerically-calculated phase of the nonlinear susceptibility, which is the same 
regardless of window width. The difference between the retrieved phase and the exact phase is shown in Figure SI 
d. As seen there are no remnants of Raman peaks; thus, this remainder is additive and represents the error term, e. 

SI.2 Comparison of Previous Derivations of the Kramers-Kronig Relation to Coher¬ 
ent Raman Spectra 

Previous use of the KK[S1, S2] utilized amplitude step-functions in the time-domain to effectively perform the Hilbert 
transform. Below we will show that these methods are, in fact, equivalent in phase retrieval and both perform a 
Hilbert transform. These derivations have been altered to conform to the notation used within this manuscript and 
a particular Fourier transform convention. Accordingly, the derivations we present may have different multiplicative 
constants, but the fundamental results are identical to the original derivations. 

SI.2.1 The “Time-Domain Kramers-Kronig Transform” (TDKK) 

Liu, Lee, and Cicerone[Sl] developed a time-domain Kramers-Kronig transform with a mathematical relaxation of 
causality. Chiefly, the Fourier-transformed (time-domain) CARS spectrum is cut at t < 0 and the time-domain NRB 
spectra is conversely cut at t > 0. Explicitly: 

(j){uj) = - 2Im|v>{ln[/cAfls(a;),/AfflB(u;)]} - ^"^^CARsUj) | ^ 

*To whom correspondances should be address 
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Figure SI: Window edge effect produces linear, additive error, (a) Simulated XR xnr- The simulation window 
is broad enough to always capture the Raman peaks, (b) Phase retrieved via the windowed Hilbert transform over 
varying spectral windows, showing good agreement in peak amplitudes and shapes but with varying baselines, (c) 
Exact phase of XR + Xnr- (d) Difference of the retrieved (b) and exact (c) phases, showing only a smooth baseline. 
This indicates that the window edge effect is (or approximately) additive. 


where ‘Im’ selects the imaginary component, and "0 is an operator defined as: 

ip{\n.[IcARsi^), iNRsi^^)]} = -F {u(t)J'”^{ln/cAiJs(w)} + u(—t)J^“^{ln/Ari?B(w)}} . (S4) 
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In eq 4, u{t) is a step-function defined as 1 for t > 0. The operator ip selects the Fourier transform of the CARS 
signal for t > 0 and the Fourier transform of the NRB signal for f < 0. Combining eqs 3 and 4 and noting that 
F{u{t)} = ^jTri2\Pj(i-KiS) -I- 5(w)], where V is the Cauchy principle value: 
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where * is the convolution operation. Using the definition of the Hilbert transform for an arbitrary function f{x): 

(S6) 


Mfix)} = — 
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fix') , V 
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X — x' 7TX 


and combining with eq 5: 

(p{uj) ='H{lnIcARsiix)} - 'HilnlNRsi^x’)} = "H |ln | = 

thus, the retrieved phase is identical with that presented in the main text of the manuscript. 


SI.2.2 The “Phase-Corrected Kramers-Kronig” (PCKK) 


Masia, et al.[S2] presented a phase retrieval method in which prior to the Kramers-Kronig transform the CARS 
signal is normalized by the NRB reference spectrum. Afterwards a step-function is applied in the time-domain and 
a Fourier-transform applied: 


(p{u}) 
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iNRsil^) / / 
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Applying the definition of the Hilbert transform, eq 6, to eq 8: 


1 IcARsfxi) _ A'. [, IcARsi‘x) \ _A-. f 1 

TTOJ Ijvrb(oj) I InrbH I IIxwbMIJ 


(S8) 


(S9) 


From a phase-retrieval point-of-view, the results between the TDKK and the PCKK are identical (eqs 7 and 9) and 
with the Hilbert transform derivation presented within the main text. In application, there is a difference between 
the PCKK and the TDKK: amplitude normalization by the NRB. In the TDKK, the Raman-like spectrum is: 

iTDKKiu}) = a/TcarsM sin(/) = |Csi(w)||x^^^(w)| sin(/), (SIO) 


where Cst is the effective stimulation profile presented within the main text of this manuscript. The TDKK 
manuscript neglected the shape of the excitation sources; thus, the retrieved spectrum was directly proportion 
the spontaneous Raman spectrum. The PCKK manuscript accounted for these sources as: 


Ipckk)^) 


' IcARsj^) 

Inrb{x>) 


sin^ 


Ixjvb(w)| 


(Sll) 


thus, the stimulation profile is removed but the output spectrum is now scaled with respect to the NRB. 

Under an ideal circumstance, the stimulation profile would be directly measurable, removed, and the retrieved 
spectrum would follow that of the TDKK. In practice, however, this is not trivial. Normalization by the NRB 
signal, as presented in the PCKK and this manuscript, removes the stimulation profile and other static, spectral 
perturbations such as the optical filter passband oscillations. Thus, this practice is a near-necessity for creating 
spectra that are directly comparable from system-to-system. As addressed within the Discussion section of the main 
manuscript, this may be alleviated with future development of reference materials or other techniques of quantifying 
the optical system spectral profile. 
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S2 Materials and Methods 

S2.1 Processing Methodology and Performance 

The typical steps taken in order to extract the Raman features from BCARS spectra or images are outlined in Figure 
S2. Some of these steps are common, others less so. As such, we will detail some of these steps in further detail 
below. 



Figure S2: Flow chart detailing the method of extracting quantitatively reliable Raman-like spectra from BCARS. 
Results of certain operations have callouts to their corresponding equations within the main text. 
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S2.1.1 Residual Background Subtraction 


Typically, 100 to 1000 dark spectra are collected with no sample illumination, averaged, and subtracted from the 
BCARS and reference NRB spectra. Due to detector operating conditions or stray light, the exact level of dark signal 
may vary (slightly) from pixel-to-pixel or image-to-image. To remove this residual dark signal, we take advantage 
of the several hundred spectral pixels that ideally do not receive anti-Stokes photons. As shown in Figure S3 a, 
^350 spectral pixels (^900 cm“^) are not illuminated, from which we can evaluate the mean (see Figure S3 b) and 
subtract. Although not explicitly discussed within this manuscript, deviations of the dark level will induce amplitude 
and phase errors within the extracted Raman-like spectra that are not straight-forward to remove. 



Figure S3: Dark signal and residual subtraction, (a) BCARS and dark spectra. The acquired spectral window 
extends beyond the boundaries set by the optical filter set. (b) After dark signal subtraction, stray light or detector 
conditions may cause the baseline level to drift. This residual may be calculated from [ideally] unilluminated portions 
of the spectrometer and subtracted on a pixel-by-pixel basis. 


S2.1.2 Denoising via the Anscombe Transformation and Singular Value Decomposition (SVD) 

SVD is a matrix factorization technique with many uses, including noise reduction, and has been applied previously 
to BCARS hyperspectral imagery[S3, S2, S4, S5, S6, S7]. BCARS hyperspectral imagery is unfolded into a two- 
dimensional matrix. A, with rows representing the spectral axis and columns spatial content. The SVD algorithm 
factorizes this matrix into three components: 


A = USV*, 


(S12) 
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where in this context, U contains the spectral bases (orthonormal eigenfunctions), S is a diagonal matrix containing 
the “singular values” (SV) in descending order (descending average contribution), V describes the spatial distribution 
of the bases in U, and is the conjugate transpose. A simplified and typical method of de-noising is to analyze the 
normalized intensity of the singular values (S), select a cut-off (C), set all higher diagonal elements S to 0, and to 
construct a denoised hyperspectral data matrix, Adenoise, as: 

^denoise = US{l:C'}VA (S13) 

This method assumes that the signal and noise are separable and that (a) the signal is entirely enveloped in the 
lowest SVs and (b) the signal will be contained in consecutive SVs below a certain cut-off. Whether these conditions 
are met is determined by the signal-to-noise ratio of the signal and the noise distribution statistics. Of significant 
consequence, is that in its most general form, SVD (and the related principle component analysis [PCA]) assume 
the noise is additive and follows a normal distribution. In BOARS, and many other techniques, the noise of often 
of a mixed nature: containing (approximately) additive white Gaussian noise (AWGN) and Poisson noise. This 
is especially important in view of the recently-developed BOARS system, which generates spectra covering a large 
intensity range. This leads to mixed-noise with the mixing ratio varying dramatically across each individual spectrum. 
There is, therefore, a necessity to “whiten” the noise as to be approximately constant (statistically) across each 
individual spectrum. 

One strategy is to perform “variance-stabilization” using an Anscombe transformation's, S9]. Figure S4 a shows 
the results of 1000 simulations of a broad spectral peak containing AWGN (standard deviation, a, of 10) or mixed 
noise. The level of the signal and the AWGN approximates the intensities found with actual BOARS experiments. 
Within spectral regions with weak intensity, the AWGN is the dominant noise source. At higher intensities, however, 
the Poisson noise is up to 10 x larger, as shown by the plots detailing standard deviation in Figure S4 b. Applying 
the generalized Anscombe transformation to the 1000 simulated spectra and re-evaluating the standard deviation, 
as shown in Figure S4 c, there is no noticeable variation of the noise across the spectra. 

Figure S5 a shows the normalized contribution of each spectral basis function (i.e., S'm.m/during 
pre-processing of the murine pancreas tissue presented within the main text. As previously discussed, the typical 
approach to denoising is to set cutoff value based on singular value contribution (such as in Figure S5 a) or via 
qualitative analysis of the spatial and/or spectral components (V and U, respectively). Figures SS5 b-g shows the 
spatial distribution of selected SVs (V-components). The spatial distributions of the lowest SVs clearly show features 
of the tissue, but by SV 13 or 14, it becomes less clear. Although there is no obvious spatial content of SV 14, the 
spectral basis function (Figure S5 h) shows Raman features. Additionally, the noise-like features below 1000 cm“^ 
are due to the Poisson noise which is largest at the lowest wavenumbers in this particular BOARS system. Figure 
S5 i shows the 100th spectral basis function and, as expected, the Poisson noise has moved to lower wavenumber 
regions as this contribution is smaller. Additionally, one can see remaining Raman features near ~2900 cm“^. This 
supports the previous supposition that with the mixed-noise of BOARS spectra, SVD does not properly separate 
signal and noise, which also results in hundreds of SVs necessary to capture the signal accurately. 

Figure S6 a-f shows the spatial distribution of SVs when using the Anscombe transformation on the BOARS 
spectra prior to SVD. In this case, SV 24 shows small spatial content and SV 25 none obvious. The spectral bases for 
these SVs (Figure S6 g,h) demonstrates that SV 24, which has spatial features, also contains spectral features, and 
SV 25 which has no spatial content also shows no obvious spectral content. Thus, there appears to be correlation 
between the spatial and spectral features. Additionally in both of these cases, there is no obvious noise variation 
across the spectrum (in contrast to Figure S5 h,i). Figure S6 i,j shows the spectral bases for the 50th and 100th 
SV. The non-variance stabilized SVD showed clear spectral features at the 100th SV, but using the Anscombe 
transformation, even the 50th SV appears to only contain noise. 

Although this section qualitatively indicates that the Anscombe transformation assists the SVD in properly 
separating noise and signal, there is no clear indication that there exists a defined boundary between SVs containing 
predominantly signal and those predominantly noise. Additionally, the task of identifying SVs containing signal 
often involves human intervention and qualitative analysis. In the next section, we present a developed method of 
automating spectral and spatial feature analysis. As will be demonstrated, improper selection of SVs may result in 
distorted spectra, with improper Raman-peak amplitudes and even extra (or missing) Raman peaks. 

S2.1.3 Automated Singular Value Selection 

Selecting the proper combination of SVs to reconstruct the original Raman content is of vital importance. Although 
the use of few SVs generates clear, attractive imagery, the underlying data may be distorted or completely erroneous. 
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Figure S4: Simulated variance stabilization via the Anscombe transformation, (a) Distribution of simulated signal 
over 1000 simulations with additive white Gaussian (AWGN) and mixed noise, (b) Standard deviation showing the 
intensity dependence of mixed noise, (c) Distribution of Anscombe transformed mixed noise signals, (d) Standard 
deviation showing normally distributed noise. 


Figure S7 a shows three spectra from within the internal elastic lamina of the artery (primarily composed of elastin) 
imaged within the main text of the manuscript. Each spectrum is the mean from the same 10 pixels but with varying 
number of SVs used in de-noising (or without SVD). The use of three SVs, shows significant spectral distortion. The 
spectral profile within the CH-stretch region would seem to indicate that the elastic lamina contains a significant 
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Figure S5: Effect of mixed noise on SVD. (a) Content fraction of each basis spectrum with and without the Anscombe 
transformation, (b-g) Spatial distributions of selected SVs. (h) The 14th basis function shows spectral features and 
spectrally narrow noise (due to Poisson noise) though (g) does not indicate a contribution, (i) The 100th basis 
function still shows some spectral features and confined noise. With mixed noise characteristics of the measured 
spectra, SVD does not well-separate signal and noise components. 


lipid content (based on the pronounced 2850 cm“^). Excluding baseline drift, the use of the first 100 SVs presents a 
significantly improved spectrum but with reduced efficacy of noise suppression (especially for a single pixel, as shown 
in Figure S7 b). This indicates that there is a clear need to select enough SVs to maintain spectral integrity and few 
enough to provide significant noise reduction. The developed method uses both spectral and spatial features and the 
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Figure S6: SVD of Anscombe transformed spectra, (a-f) Spatial distribution of selected SVs (Ansc SV). (g) The 24th 
spectral basis function shows clear features as would be indicated by spatial content in (e). (h) The 25th SV appears 
to only contain noise (as does the corresponding (f)). By the 50th (i) and 100th (j) SVs, there is no noticeable 
remnants of signal. Without the Anscombe transformation, even the 100th SV (Figure S5 i) still contained Raman 
peak information. 


Fourier-domain noise statistics. 

Figure S8 a,b demonstrates the concept of spatial domain SV selection. The absolute value of the spatial 
components (mean-subtracted) for each SV (|V|) are transformed via a two-dimensional Fourier transform. A 
rectangular boundary is applied to separate the low-frequency components (assumed to contain primarily signal) 
and high-frequency components (assumed to be primarily noise). The “spatial signal ratio” is defined as the sum 
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Figure S7: Raman-like spectral distortions from improper singular value selection, (a) Mean spectra of 10 pixels 
within the elastic lamina of the murine pancreas tissue. Using only the first 3 SVs results in significant errors. Using 
the first 100 SVs shows significant improvement, (b) Single-pixel spectra from within the elastic lamina showing 
similar distortions with too few SVs incorporated. 
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of pixel intensities (complex modulus of the Fourier transform components) within the boundary divided by those 
outside the boundary (scaled by the number of pixels within and outside of the boundary). Figure S8 c plots the 
spatial signal ratio as a function of SV number. As expected the lowest SVs show the largest ratio. To develop 
a cut-off value, we assume that the highest SVs (>700) contain only noise. From the spatial signal ratio of these 
highest SVs, we calculate the standard deviation, cr. The cut-off value, for which all SVs below will be neglected, is 
heuristically determined as a multiple of cr. Based on visual inspection of V, we select a multiplier of 3.5. This cut-off 
value may not be optimal, but it does appear, in our experience, fairly robust across images of different samples; 
although, we perform manual inspection of selected SVs prior to further progress in the pre-processing workflow. 
Figure S8 c shows this cut-off value in red. Between SV 15 and 20 (see Figure S8 c, inset), the spatial ratio drops 
below the cut-off (SV 18) then returns above (SV 19). Images of the spatial distribution of SV 18 and 19 (Figure S8 
c) demonstrate that our method correctly categorized SV 18 as containing no obvious spatial components but SV 
19 does. This demonstrates that there may be SVs with primarily noise components intermixed between SVs with 
primarily signal components. 

As illustrated in Figure S5 g,h, the lack of spatial components may or may not indicate a lack of spectral 
components. Thus, we have also developed an automated spectral basis set analysis technique. As pictographically 
described in Figure S9 a,b, the spectral component selection tool analyzes the statistics of the U components within 
the Fourier domain (time-domain). Again, a ratio is determined between the sum of the signal contribution within 
a low-frequency window and outside this window (scaled to number of pixels). The half-width of this window is 
set to correspond to the temporal duration of our probe source (~ 3.4 ps). Again, we assume that the highest SV 
components contain primarily noise and calculate the standard deviation. The spectral components, however, show 
a mean-line drift, which we fit with a low-order polynomial (in this case, third-order). The cut-off is a multiple 
(4, in this case) of the standard deviations from this mean. Figure S9 c (inset) demonstrates qualitatively correct 
identification of SV 26 with no noticeable Raman components and SV 27 with signal features. 

Table SI shows the time of each pre-processing step in automatically processing the murine pancreas tissue 
presented within the main text. The computer was a Dell Optiplex 9010 with quad-core Intel 17-3770 CPU at 
3.4 GHz, with 16 GB of memory, and running MATLAB R2013a. The total processing time was approximately 
28 minutes. At 95% of the computation time, the most intensive process was the automated detrending using 
asymmetric least squares (ALS). This method, although robust and reuqiring no user intervention, requires multiple 
matrix inversions, which do not facilitate processing several spectra in parallel. It should be noted that the ALS 
algorithm, as developed, uses the CHOLMOD implementation of fast sparse matrix inversion for optimal performance. 
Future work will investigate alternative means of automated detrending. Excluding the ALS step, the total processing 
time was approximately 90 seconds (< 1 ms/spectrum). 


Sub-Process 

Total Computation Time (s) 

Time/Pixel (ms) 

Dark Subtraction 

0.74 

0.01 

Residual Baseline Subtraction 

1.14 

0.01 

Anscombe Transform 

0.48 

0.01 

SVD 

41.71 

0.46 

Automated SV Selection 

3.67 

0.04 

Inverse Anscombe Transform 

5.18 

0.06 

Phase Retrieval (KK) 

17.66 

0.20 

Phase Detrending (ALS) 

1564.91 

17.44 

Spectral Scaling 

19.69 

0.22 

Total 

1655.18 

18.45 


Table 1: Computation time of automated pre-processing of hyperspectral BCARS image of murine pancreas. The 
total processing time was ~28 minutes. 


For comparison with the demonstrations above. Figures SSIO and SSll show the automated spatial and spectral 
analysis results for non-Anscombe transformed BCARS data. 
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Figure S8: Automated selection of singular values based on spatial content. (a,b) The spatial content is converted 
into the Fourier-domain from which “signal” and ’’noise” regions are defined, (c) The spatial signal fraction is 
proportional to the total intensity ratio of the “signal” and ’’noise” regions. A cut-off is selected based on a user- 
selected multiplier of the standard deviation of the spatial signal ratio at the highest SVs, which are assumed to 
predominantly contain noise. 


S3 Results 

S3.1 The Maximum Entropy Method (MEM) and Phase Error Correction Applied 
to Simnlated Spectra 

The MEM performs phase-retrieval based on information theory ground[S10]. Previously, the MEM and KK methods 
were demonstrated to be “functionally” equivalent[S11]. To that end, we performed phase retrieval via the MEM on 
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Figure S9: Automated selection of SVs based on spectral content. (a,b) The spectral basis functions are converted into 
the time-domain from which “signal” and ’’noise” regions are defined, (c) The spectral signal fraction is proportional 
to the ratio of the “signal” and ’’noise” regions. A cut-off is selected based on a user-selected multiplier of the 
standard deviation of the spectral signal ratio at the highest SVs, which are assumed to predominantly contain noise. 


simulated BCARS spectra containing two peaks under ideal (known NRB) conditions and with a surrogate reference 
NRB. As shown in Figure S12 a, the retrieved Raman spectra differ between the MEM and KK methods, with the 
MEM showing a slightly larger baseline drift. Upon phase detrending and scaling, the Raman-like spectra extracted 
using reference spectra agree with those with the known NRB. Differences in peak amplitudes between the MEM 
and KK (<5%) are due to the underlying algorithmic differences between the KK and MEM methods. Comparison 
of these methods is presented in further detail in Ref. Sll. 
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Figure SIO: Automated selection of singular values based on spatial content without use of the Anscombe transfor¬ 
mation. In comparison to the results in Figure S8, without the use of the Anscombe transformation, SVs showing 
clear spatial contributions are spread over several hundreds SVs, such as 117 and 177 (inset). This automated method 
improves SV-selection time as it would be laborious for a user to manually inspect hundreds of SVs. Additionally, 
the automated method removes many of the SVs that contain primarily noise; although, they are at lower SVs. 


53.2 BCARS Spectra of Glycerol 

As described within the main text, glycerol spectra were collected on two separate BCARS platforms and pre- 
processed with different reference NRB spectra. Figure S13 a shows BCARS spectra of glycerol collected on a 
recently-developed platform[S3] (“System 1”) and an older system[S12, S6] (“System 2”). There is a clear difference 
in the system responses with the newer system showing marked enhancement at the lowest wavenumbers. Figure 
S13 b shows the BCARS spectra of 3 different reference materials (and coverslip glass acquired on System 2). These 
spectra show differences in amplitude and shape, which will be discussed in detail in the next subsection. Figure 
S13 c shows the retrieved Raman spectra using the KK relation (no error correction performed). Inspection of the 
Raman peaks (such as the CH-stretch peaks) indicates there is substantial differences in retrieved amplitude and 
phase. 

53.3 Time-Window Self-Referencing to Directly Capture the Approximate NRB of 
the Sample and the Raman Peaks of Reference Materials 

Time-window self-referencing (TWSR) is a method we have implemented to capture a CARS spectrum that predom¬ 
inantly contains the NRB with reduced contributions of the Raman vibrational components. This techniques can be 
used to retrieve an approximate NRB spectrum on a pixel-by-pixel basis, requiring a second image to be acquired, 
but this falls outside of the scope of this manuscript. Rather we will use TWSR to enable extraction of the Raman 
signature of typical NRB reference materials to highlight that these materials do, in fact, contain Raman peaks that 
will affect the final spectral retrieval. 

As described within the main text, the “effective” nonlinear susceptibility, y, is related to the nonlinear suscep¬ 
tibility as: 

xi^) = xi^)*Epriuj), (S14) 

where * is the convolution operation and Epr is the electric field of the probe source. In essence the spectral resolution 
of the CARS spectrum is determined by the spectral narrowness of the probe source [SI 3]: the spectrally narrower 
the probe source, the higher the resolution of the CARS spectrum. From a time-domain perspective, this can be 
viewed as how much temporal information is acquired. One can describe the effective time response of the nonlinear 


14 




















Without Anscombe Transformation 



Singular Value Number 

Figure Sll: Automated selection of singular values based on spectral content without use of the Anscombe transfor¬ 
mation. In comparison to the results in Figure S9, without the use of the Anscombe transformation, SVs showing 
clear spectral contributions are spread over several hundreds SVs, such as 149 and 229 (inset). Additionally, without 
the Anscombe transformation, there is an obvious mixing between signal contributions and those due to Poisson 
noise. 


susceptibility, R{t), as: 

m = {Epr{uj)h (S15) 

^ V 

R{t) -Epr(t) 

where is the inverse Fourier transform, R{t) is the time response of the material nonlinear susceptibility, and 
Epr{t) is the temporal field profile of the probe source. Figure S14 a shows the time response of a simulated nonlinear 
susceptibility (resonant and nonresonant contributions) containing two Raman peaks and a broad (25,000 cm“^ full- 
width half-max) nonresonant background (see Figure SI a for plots of the individual terms in the frequency-domain). 
Notice that the nonresonant term only contributes a signal for a brief duration (femtoseconds) but the Raman 
contributions decays over several picoseconds. Under normal spectroscopic collection, as described in Figure S14 
b, the temporal overlap of the material stimulation (via pump and Stokes sources) and the probe source is set to 
maximize the total energy collection (i.e., from time 0, on); thus, maximizing the temporal duration captured and 
the spectral resolution. In the TWSR technique, the probe source is temporally offset as to only capture the earliest 
moments of signal creation (Figure S14 c); thus, acting as a femtosecond probe. The signal generation from the 
nonresonant component of the nonlinear susceptibility is the same as under normal operating conditions, but little 
of the Raman decay is acquired. The “early time” (ET) spectrum will serve as an approximate NRB measurement. 

As an experimental demonstration of the effect of Raman peaks within common reference materials. Figure S15 a 
shows the retrieved Raman-like spectra of glycerol within the lower wavenumber region using reference NRB (normal 
temporal probe settings). The spectra show distinct similarity, but there are certain regions of deviation at ^560 
cm“^, ^920 cm“^, and ^1100 cm“^. By using the spectra gathered from these reference materials and the TWSR 
technique to capture an approximate NRB spectrum, one can retrieve Raman-like spectra from the reference materials 
themselves as shown in Figure S15 b. Comparing Figure S15 (a,b), it is apparent that the Raman peaks of coverslip 
and glass slide correspond to the same spectral regions of deviation. Figure S15 c shows the retrieved Raman-like 
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Figure S12: Phase detrending and scaling applied to the maximum entropy method (MEM), (a) Comparison of 
extraction of Raman signatures using the KK and MEM methods with a known NRB (“ideal”) and with a reference 
NRB (“Ref’). Both methods produce distorted spectra due to phase and amplitude errors when using a reference 
NRB. (b) Applying the previous described phase error corrections technique produces respective spectra for the 
MEM and KK that are identical whether using the actual NRB or a surrogate reference. Differences between the 
MEM and KK are due to the underlying numerical algorithms and not explicitly to due to phase or amplitude errors. 


spectra of glycerol using TWSR with the reference materials, reducing their Raman spectral contribution. The use 
of TWSR in this manner reduces the deviation of spectral components by a factor of ~2. 

S3.4 Tissue Imaging and Histogram Analysis 

The pseudocolor imagery and analysis was performed to highlight general protein content, smooth muscle (predom¬ 
inantly actin/myosin), and DNA/RNA. For protein 2937 cm“^ - 2882 cm“^ was used, which effectively suppresses 
the lipid content. For smooth muscle, a peak at 1339 cm“^ was used. As this particular peak is a portion of a 
shoulder that is spectrally broad, containing several neighboring peaks, a linear interpolant was calculated between 
1288 cm“^ and 1360 cm“^, subtracted, and the peak amplitude at 1339 cm“^ was calculated. For DNA/RNA, the 
peak at 785 cm“^ was used. To mitigate the effect of residual baseline, this peak amplitude was calculate relative to 
a linear interpolant calculated between the neighboring troughs at 763 cm“^ and 820 cm“^. 

As described in the main text, phase detrending and scaling can significantly reduce the effect of using reference 
NRB spectra. Figure S16 a-c shows histograms of the calculated peak amplitudes relating to protein, smooth muscle, 
and DNA/RNA as a function of reference NRB material corrected with traditional amplitude detrending. For the 
case of general protein (Figure S16 a), the use of glass coverslip as the surrogate NRB retrieves peak amplitudes 
with a larger mean and broader distribution than with water. The converse is true for smooth muscle (actin/myosin) 
for which retrieval increases the amplitude distribution to higher values (Figure S16 a). For DNA/RNA, the re¬ 
trieved amplitudes are relatively similar. Figure S16 d-f demonstrates that the use of phase detrending and scaling 
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Figure S13: Extraction of Raman-like spectra of glycerol using difference NRB references and difference BCARS 
platforms, (a) Raw BCARS spectra collected ona recently-developed BCARS platform (“System 1”) and a more 
traditional implementation (“System 2”), showing significantly difference system responses, (b) Raw BCARS spectra 
of 3 common NRB reference materials: glass coverslip, glass microscope slide, and water. Additionally, a raw spectrum 
from glass coverslip was collected on System 2. (c) Raman-like spectra extracted from (a) using spectra in (b) without 
amplitude or phase error correction, showing different shapes and amplitudes for corresponding peaks. 
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Figure S14: Time-window self-referencing (TWSR) to retrieve the approximate NRB of the sample, (a) Temporal 
response of the material calculated from the fast Fourier transform (FFT) of the nonlinear susceptibility, (b) Under 
normal operation, the probe source delay (multiplicative in the time-domain) is set to acquire the entire temporal 
evolution of the Raman/CARS signal, (c) One can acquire a spectrum that is predominantly NRB by setting the 
probe delay to an early time point. There still exists a Raman contribution, but it is weaker and the peaks are 
significantly broadened owing the effective, limited spectral resolution. 


significantly improves the similarity of retrieved amplitudes regardless of reference material used. 


Disclaimer 

Any mention of commercial products or services is for experimental clarity and does not signify an endorsement or 
recommendation by the National Institute of Standards and Technology. The authors declare no competing financial 
interests. 
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Figure S15: Using temporal dynamics of the CARS generation process to improve the utility of NRB reference 
measurements, (a) The retrieved Raman-like spectra of glycerol (with the use of reference NRB measurements 
and corrected via phase detrending and scaling) showing close agreement but some deviations, (b) Raman-like 
spectra extracted from common reference NRB materials using TWSR to collect an approximate NRB spectrum of 
each material. The peaks in (b) correlate with the differences between spectra in (a), (c) Extracted Raman-like 
glycerol spectra using NRB reference spectra that were collected under early-time conditions to reduce the Raman 
contributions. Comparing (a) and (c), the spectral differences are reduced in amplitude by ~2x. 
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Figure S16: Histogram analysis of spectral peaks amplitudes using traditional amplitude detrending and the new 
phase detrending methods, (a-c) Pixel-by-pixel comparison of Raman-like peak intensities using different NRB surro¬ 
gate materials and traditional amplitude detrending, (d-f) Pixel-by-pixel comparison of Raman-like peak intensities 
using different NRB surrogate materials and the developed phase detrending and scaling method. Note: histogram 
bins are always the same width- bar widths adjusted for visual clarity. 
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